The program is written in R.
View Code
function(input_type, nsim, n, lambda1, hratio, hratio_alt, accrual, follow,
event_frac, int_frac, int_follow, alpha, alpha_interim, frac_control) {
eventv1 = eventv2 = int.sev = int.eventv1 = int.eventv2 = int.betav = betav = testv = rep(NA, nsim) # storage
lambda2 = lambda1 / hratio # hazard rate for the experimental
int_accrual = min(int_frac, 1) * accrual # int_frac must be smaller or equal to 1
sample.size = n
n1a = round(frac_control * n)
n2a = sample.size - n1a
# If provided, convert event fraction to accrual and follow-up time
if (input_type == "Total") {
follow.v = seq(from = 0, to = follow, by = follow / 100)
accrual.v = seq(from = 0, to = accrual, by = accrual / 100)
prob1.arm1 = (accrual.v - 1 / lambda1 + (1 / lambda1) * exp(-lambda1 * accrual.v)) / accrual #arm 1 prob event
prob2.arm1 = 1 - exp(-lambda1 * follow.v) + (exp(-lambda1 * (accrual + follow.v)) / lambda1 +
(accrual - 1 / lambda1) * exp(-lambda1 * follow.v)) / accrual
prob1.arm2 = (accrual.v - 1 / lambda2 + (1 / lambda2) * exp(-lambda2 * accrual.v)) / accrual #arm 2 prob event
prob2.arm2 = 1 - exp(-lambda2 * follow.v) + (exp(-lambda2 * (accrual + follow.v)) / lambda2 +
(accrual - 1 / lambda2) * exp(-lambda2 * follow.v)) / accrual
probv.arm1 = c(prob1.arm1, prob2.arm1)
probv.arm2 = c(prob1.arm2, prob2.arm2)
probv = frac_control * probv.arm1 + (1 - frac_control) * probv.arm2 # sum events accross arms weighted by arm sample size
times = c(accrual.v, follow.v)
accstat = c(rep(1, 101), rep(2, 101))
ind1 = probv >= event_frac
indpt = min((1:202)[ind1])
if (accstat[indpt] == 1) {
int_accrual = times[indpt]
int_follow = 0
} else if (accstat[indpt] == 2) {
int_accrual = accrual
int_follow = times[indpt]
}
int_frac = min(int_accrual / accrual, 1)
}
# This is the simple case
for (jj in 1:nsim) {
# Generate full analysis data unobservables
tt1 = -log(runif(n1a)) / lambda1
tt2 = -log(runif(n2a)) / lambda2
cc1 = follow + runif(n1a) * accrual
cc2 = follow + runif(n2a) * accrual
tt = c(tt1, tt2)
cc = c(cc1, cc2)
# Generate times corresponding to the first interim
int.n1 = round(int_frac * n1a) #no pts in arm 1 at interim
int.n2 = round(int_frac * n2a)
oo1 = rev(order(cc1)) #ordering for accrual time
int.patients1 = oo1[1:int.n1]
oo2 = rev(order(cc2))
int.patients2 = oo2[1:int.n2]
int.cc = c(cc1[int.patients1], cc2[int.patients2]) - (follow - int_follow) #create true censoring times
int.tt = c(tt1[int.patients1], tt2[int.patients2]) #create true survival times
int.cc = int.cc - min(int.cc) + .01 #??? Why wouldn't I do the same for true survival times
int.status = 1 * (int.tt < int.cc)
int.group = c(rep(0, int.n1), rep(1, int.n2))
int.times = int.status * int.tt + (1 - int.status) * int.cc
int.bb = survival::coxph(survival::Surv(int.times, int.status) ~ int.group) # So we can store all the parameters we need from here
int.betav[jj] = int.bb$coef
int.sev[jj] = summary(int.bb)$coef[3]
int.eventv1[jj] = sum(int.status[int.group == 0])
int.eventv2[jj] = sum(int.status[int.group == 1])
# final analysis
status = 1 * (tt < cc)
group = c(rep(0, n1a), rep(1, n2a))
times = status * tt + (1 - status) * cc
bb = survival::coxph(survival::Surv(times, status) ~ group)
testv[jj] = bb$score
betav[jj] = bb$coef
eventv1[jj] = sum(status[group == 0])
eventv2[jj] = sum(status[group == 1])
}
# mean number of events in control/experimental arm at time of each analysis
events.int = c(mean(int.eventv1), mean(int.eventv2))
events.final = c(mean(eventv1), mean(eventv2))
event_mat = rbind(c(events.int, sum(events.int)), c(events.final, sum(events.final)))
rownames(event_mat) = c("interim", "final")
colnames(event_mat) = c("Control", "Expt", "Total")
vv = dim(event_mat)
# alpha = .2 # Testing Null
# alpha_interim = .05 # Testing Alternative
# Calulate Futility Stopping Probablities
# Using HR = 1 Rule
stop.betav =! (int.betav < 0)
stop.HR1 = mean(stop.betav)
# Using Test with alpha_interim
tt2 = (-((-log(hratio_alt) - int.betav) / int.sev)) < qnorm(1 - alpha_interim)
stop.testv =! tt2
stop.test = mean(stop.testv)
# Thresholds Using Test
stop.test.thresh = -(log(hratio_alt) - qnorm(1 - alpha_interim) * int.sev)
# Power with No Interim Analysis
rejectv = (testv * (betav < 0)) > qchisq(1 - alpha, 1)
power.no = mean(rejectv)
# Power with Interim Analysis - using HR1 rule
power.HR1 = mean(rejectv * (!stop.betav))
# Power with Interim Analysis - using Interim Test
power.test = mean(rejectv * (!stop.testv))
# OUTPUT
output =
list(n_per_arm = n,
lambda1 = lambda1,
hratio = hratio,
accrual = accrual,
follow = follow,
int_frac = int_frac,
int_accrual = int_accrual,
int_follow = int_follow,
alpha_interim = alpha_interim,
alpha = alpha,
total_events = event_mat[vv[1], vv[2]],
event_mat = as.data.frame(event_mat),
n_interim_analysis = int.n1,
p_stop_int_hr = stop.HR1,
p_stop_int_test = stop.test,
p_stop_int_test_minus_beta = mean(stop.test.thresh),
power_no_minus_interim = power.no,
power_hr1_minus_interim = power.HR1,
power_test_minus_interim = power.test)
return(jsonlite::toJSON(output, pretty = TRUE))
}